//kirk.cpp

#include <math.h>
//#include <stdio.h>

static double max,f1,f2;
double mxf(double a,double b);
void hsu0_init(double min1,double max1,double f10,double f20){

//***************** Parameter initialization **********************************

      f1=f10;
      f2=f20;
      max=mxf(fabs(min1),fabs(max1));
return;}

//-----------------------------------------------------------------------------

double hsu0(double u,double v,double a){

//********* Phase function computation at single point ************************
      double psi,psi0=0;
      if(a<0) psi0=M_PI;
      psi=psi0+acos(fabs(a)/max)*pow(-1,int(f1*u+f2*v));
return psi;}

void hsu1_init(double min1,double max1,double f10,double f20){

//***************** Parameter initialization **********************************

      f1=f10;
      f2=f20;
      max=mxf(fabs(min1),fabs(max1));
return;}

//-----------------------------------------------------------------------------

double hsu1(double u,double v,double a){

//********* Phase function computation at single point ************************
      double psi,psi0=0;
      if(a<0) psi0=M_PI;
      psi=psi0+asin(fabs(a)/max)*pow(-1,int(f1*u+f2*v));
return psi;}

